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ABSTRACT 

We present our numerical results of two-dimensional magnetohydro dynamic 
(MHD) simulations of the collapse of rotating massive stars in light of the collap- 
sar model of gamma-ray bursts (GRBs). Pushed by recent evolution calculations 
of GRB progenitors, we focus on lower angular momentum of the central core 
than the ones taken mostly in previous studies. By performing special relativistic 
simulations including both realistic equation of state and neutrino cooling, we 
follow a long-term evolution of the slowly rotating collapsars up to ~ 10 s, ac- 
companied by the formation of jets and accretion disks. Our results show that for 
the GRB progenitors to function as collapsars, there is a critical initial angular 
momentum, below which matter is quickly swallowed to the central objects, no 
accretion disks and no MHD outflows are formed. When larger than the crite- 
ria, we find the launch of the MHD jets in the following two ways. For models 
with stronger initial magnetic fields, the magnetic pressure amplified inside the 
accretion disk can drive the MHD outflows, which makes the strong magnetic ex- 
plosions like a 'magnetic tower'. For models with weaker initial magnetic fields, 
the magnetic tower stalls first and the subsequent MHD outflows are produced 
by the magnetic twisting of the turbulent inflows of the accreting material from 
the equatorial to the polar regions. Regardless of the difference in the formation, 
the jets can attain only mildly relativistic speeds with the explosion energy less 
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than lO^^erg. To obtain stronger neutrino energy depositions in the polar funnel 
regions heated from the accretion disk, we find that smaller initial angular mo- 
mentum is favorable. This is because the gravitational compression makes the 
temperature of the disk higher. Due to high neutrino opacity inside the disk, 
we find that the luminosities of v^. and become almost comparable, which is 
advantageous for making the energy deposition rate larger. We discuss how the 
energy deposition can be as efficient as the magnetically-driven processes for en- 
ergetizing jets. Among the computed models, we suggest that the model with 
the initial angular momentum of j ~ l.Sjiso (jiso: the angular momentum of the 
last stable orbit) and with initial magnetic field strength of ~ 10^° G, provides 
a most plausible condition for making fireballs for GRBs, because such model is 
appropriate not only for producing the MHD outfiows quickly by the magnetic 
towers, but also for obtaining the stronger neutrino heating in the evacuated 
polar funnel. 

Subject headings: coUapsar: black hall, disk — supernovae: collapse, rotation — 
magnetars: pulsars, magnetic field — methods: numerical — MHD — special 
relativity — gamma rays: bursts 



Introduction 



Gamma-ray bursts (GRBs) are one of the most energetic phenomena in the universe. 
Accumulating observations discovered such as by SwifS and HETE-2% show that GRBs are 
basically categorized into two, namely short-hard and long-soft bursts (e.g., iNakarl (120071 ) 



for review). M ore surprisingl y , GRBs with some mix e d fea tures of the two types have been 
reported (e.g., iGehrels et al.l (120061 ). iGal-Yam et al.l (120061 )). The mystery of their central 
engines seems to be thickening, w hich have l o ng pu zzled astrophysicists since the accidental 
discovery in the late sixties (e.g., iMeszarod (120061 ) . for review). Speaking about the long- 
duration GRBs, a number of their host galaxies have been identified as metal poor galaxies 
(jSavaglio et al.ll2006l : IStanek et al.ll2006l . and reference therein) . The preponderance of short- 
lived massive star formation in such young galaxies, as well as the identification of SN Ib/c 
light curves in some bursts, has pr ovided strong support to identify a massive stellar collapse 
as an origin of th e long GRBs (jPaczynskil Il998l : iGalama et al.l Il998l : iHjorth et al.l 12003 



Stanek et al.ll2003l ). The duration of the long bursts may correspond to the accretion of debris 
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falling into the central black hole (BH) (jPiro et al.lll998l ). which suggests the observational 
consequence of the BH formation likewise the supernova of neutron star formation. Pushed 
by those observations, the so-called coUapsar has received quite some interest for the cen tral 
engines 



i Observations, tne so-called collapsar nas received quite some interest tor tne cen 
of the long GRBs JWooslevlll993l : IPaczvnskil EoQsl : iMacFadven fc WooslevlEoOol ) 



In the collapsar model, the central cores with significant angular momentum collapse 
into a BH. Neutrinos emitted from the accretion disk around the BH heat the matter of 
the funnel region of the disk, to launch the GRB outflows. The relativistic outflows are 
expec ted to ultimately form a fireball, which is good for explaining the observed afterglow 
(e.g., iPiraru (119991 )). In addition, it is suggested that the strong magnetic fields in the cores 



of order of lO^^G play also an active role both for driving the magn eto-driven jets and fo r 
extracting a signi f icant a mount of energy from the central engine (e.g. , I Wheeler et al.l (120001 ) ; 
Thompson et al.l (|2004| ) ; lUzdensky fc MacFadyenI (j2007al ) and see references therein) . 



To obtain clear understanding of such scenarios, it is ultimately necessary to perform 
stellar core-collapse simulations, which trace all the phases in a consistent manner starting 
from the stellar core-collapse, core-bounce, shock-stall, stellar explosion (phase 1) or BH 
formation and the formation of accretion disk (phase 2), energy deposition to the funnel 
region by neutrinos and/or magnetic fields (phase 3), to the launching of the fireballs (phase 
4). Here for convenience we call each stage as phase 1, 2, etc. The requirement for the 
numerical modeling to this end is highly computationally expensive, which necessitates the 
multidimensional magnetohydrodynamic(MHD) simulations not only with general relativity 
for handling the BH formation, but also with the multi-angle neutrino transfer for treating 
highly anisotropic neutrino radiation from the disks. Thus various approximate approaches 
have been undertaken. All the studies, which we will mention below, are complimentary in 
the sense that the different epochs are focused on, with the different initial conditions for 
the numerical modeling being taken. 

As for the phase 1, the roles of rapid rotation and magnetic fields have been elaborately 



investigated to study t he formation q 



Takiwaki et al. (2004) 



mag n etars and its impli c ations to the collapsa. r s (e.g 



Kotake et all (12004 ):ISawai et al.l ( 20051) .lObergaulineer et all (120061): 



Suwa et al. 



(2007) 



Burrows et all (120071 ): 



Takiwaki et al 



(120091 ) and collective references in 



Kotake et al.l (120061 )). After the failed or weak explosion in the postbounce phase, the ac- 
cretion to the central objects may lead to the formation of a BH (phase 2), which severa l 
general relativistic studies have focused on (IShibata et al.ll2006l : ISekiguchi fc Shibatal 120071 ). 
Treating the BH as an absorbing boundary or using the fixed metric approaches, the numer- 
ical studies of the phase 3 are concerned with the initiation of the outflows from the funnel 
region of the disk to the acceleration of the jets as a result of the neutrino heating and/or 
MHD processes till the jets become mildly relativistic. Numerical studies of the phase 4 are 
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mainly concerned with the dynamics later on, namely, the jet propagation to the breakout 
fro m the star, when th e acceleration of the jets to the hig h Lore ntz factor is expected (see. 



e.g. 



Alov et all fl2000[ ): IZhang et all (120031 ): iMizuno et all (120071 ) and references therein) 
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genitor cores have significant angular momentum 10 cm /s) with strong magnetic fields 
(> 10^^ G), magneto-driven jets can be launched strong enough to expel the matter along the 
rotational axis within several seconds after the onset of collapse. The combination of such 
rapid rotation and strong magnetic fields is recently considered to be possible for rapidly 
rotating metal-poor sta rs, which experience the so-called chemically h omogeneous evolution 



in the main sequence (IWoosley fc Hegerl l2006l : lYoon fc Langerl l2005l ). It should be noted 



that the angular momentum of those GRB progenitors (~ 10 cm /s), albeit not a final 
answer due to much uncertainty in the ste llar mass loss, angular rn o mentum transport, an d 



magnetorotational instability (MRI) (e.g.. IVink fc de Koterl (|2005| ). iDetmers et al.l (l2008l )) 



is relatively smaller than those assumed in most of t he previous co l lapsar simulations. Em- 
ploying such a GRB progenitor, it was reported by iDessart et al.l (120081 ) based on the 2D 
radiation MHD core-collapse simulations that too much angular momentum is not favorable 
for coUapsars because the MHD explosions in the immediate postbounce phase are so strong 
that the BH formations are circumvented. 

These situations motivate us to focus on slower rotation of the central core in the con- 
text of collapsar models. As for the initial magnetic fields, we choose to explore relatively 
smaller fields (< 10^° G), which has been less investigated so far. Paying particular attention 
to the smaller angular momentum, it takes much longer time to amplify the magnetic fields 
large enough to launch the MHD jets than previously estimated. By performing special 
relativistic MHD simulations, which enable us to follow a long-term evolution over ~ 10 
s, we aim to clarify how the properties and the mechanism of the MHD jets could change 
with the initial angular momentum and the initial magnetic fields. It is noted here that 
the long-term simulation could be import ant for understanding t he X-ray flares recently dis- 
covered in a number of long GRBs (e.g., iProga fc ZhangI (120061 )). As for the microphysics, 
we include both realistic equation of state ( EOS) and neutrino coo li ng, which have been o f- 
ten neglected or oversimplified (see, however, iFujimoto et al.l (120061 ): iNagataki et al.l (120071 )). 
By doing so, we estimate the neutrino luminosities emitted from the accretion disks, and 
clarify what conditions are pivotal to make the energy depositions via neutrino pair annhi- 
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lation as efficient as the magnetically-driven processes for energetizing jets. The range of 
specific angular momentum required for the progenitors of coUapsars was predicted to be 
3 < ?"ifi(= 7/lO^^cm^ s~^) < 20 by a pioneering collapsar simulation but without MHD 



( iMacFadyen fc Woosleyi Il999l ). This is because if angular momentum is too small, mass 
element cannot stay at the last stable orbit, while if too large, mass element cannot fall 
into compact objects, and thus cannot form disk, suppressing the sufficient energy release 
for GRBs. Based on our results, we hope to understand how this criterion would change if 
MHD effects are taken into account. 

We summarize the numerical methods in ^ §3 is devoted to the initial models. The 
main results are described in §11 We summarize our results and discuss their implications 
in §3 Details of the treatments of neutrino cooling/heating in the framework of special 
relativity are given in the appendices. 



2. Numerical methods 

The results presente d in this paper are ca lculated by the MHD code in special rela- 



tivity (SR) developed by iTakiwaki et al.l (120091 ). In the following, we briefly mention the 



importance of SR for collapsar simulations and summarize shortly the numerical schemes. 

SR effects are indispensable not only to follow the propagation of GRB jets in high 
Lorentz factors, but also to follow the dynamics of infalling material, because their free- 
fall velocities and rotational velocities become close to the speed of light near the central 
compact objects. Moreover, the velocity of the Alfven waves during the jet propagation can 
be estimated as, 

VA = ~ 10^°cm/s^ ' 1 

v/p/(105g/cm3)' ^ 

where p and B are the typical density and the magnetic field near along the rotational axis. 
It can be readily inferred that the Alfven velocity can exceed the speed of light unphysically 
in the Newtonian simulation, especially for the regions where the density becomes low and 
the magnetic fields become strong. Such a situation is ubiquitous for collapsar simulations. 
Even if the propagation speeds of the jets are only mildly relativistic, we have learned that 
(at least) SR treatments are quite important for keeping the stable numerical calculations 
in good accuracy over the longer-term evolution. 



The SRMHD par t of the code is based on the formalism of iDe Villiers et al.l (120031 ) (see 



Takiwaki et al. J2009h for details). To formalize the basic equations, we need two frames; 



the laboratory frame (we call shortly as lab frame), which is the center of mass system of 
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the star, and the rest frame of the relativistic fluid. These two frames are related to each 
other with the usual Lorentz transformation. Before going to the basic equations, we write 
down the definition of the primary code variables. The state of the relativistic fluid element 
in the rest frame is described by its density, p; specific energy, e; velocity, f and pressure, p. 
Magnetic fields in the lab flame is described by the 4-vector -y/ivrfe'^ = *F'^'^Uu, where *F^'^ 
is the dual tensor of the electro-magnetic field strength tensor and Ui, is the 4-velocity. 



The basic equations of the SRMHD code can be described as follows: 

: 



9D 1 ^ , 
dt V7 



djSj - b% 
dt 



-p- 



1 



dW 



dt 



V7 



(2) 
(3) 



-t^ph{Wvkf-ib,f)da'' 



d{Wb' - Wb^v'^ 
dt 



+ dj{Wv^U -Wv'V) = 



(4) 
(5) 



fc^poi 




(6) 



where W 



D = pW, E = eW and Si = phW^Vi are the Lorentz boost factor. 



auxiliary variables correspond to density, energy, and momentum, respectively. All of them 
are defined in the lab frame. Equations. ( l2|3|4l) represent the mass, energy, and momentum 
conservations. InEq.(jll) it is noted that the relativistic enthalpy, h= {1+e/ p+p/p+\bf / p) 
includes magnetic energy. Eq.([5]) is the induction equation for the magnetic fields. In 
solving the equation, the meth od of characteri s tics is implemented to propagate accurately 
all modes of MHD waves (see iTakiwaki et al.l (120091 ) for details). 5* is the magnetic field 
in the rest frame, which is related to the one in the lab frame as = WU — Wbtb"^. Here 
bt is a time component of the 4-vector of 6^. Eq.(l6]) is the Poisson equation for the (self- 
) gravitational potential of $ pni- We eraploy the realistic equation of state based on the 



relativistic mean field theory fjShen et al.lll998l ). For lower density regime (p < 10^'^g cm~^), 
where no data is available in the EOS table with the Shen EOS, we use another EOS, which 
includes contributions from an ideal gas of nuclei, radiation, and electrons and positrons with 
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arbitrary degrees of degeneracy flBlinnikov et al.lll996l ). We carefully connect two EOS at 



p = 10 g cm for physical quantities to vary continuous in density at a given temperature. 

Ci, in Eg .(3) is the n e utrino luminosity evaluated wi th a multi-flavo r leak age scheme 
jEpstein fc Pethickl Jl98lh : iRosswog &: Liebendorfeil JiooJ; iKotake et aP JSoJ), in which 
Ue, i^e, and the heavy-lepton neutrinos, u^, z/^, Ur, (collectively referred to as ux) are 
taken into account. Electron capture on p roton and fre e nucle i, positron capture on ne utron, 
photo-pair, plasma process are included (IFuller et al.l (Il985l ). iTakahashi et al.l (119781 ) . Itoh 
et al. 1989, 1990). Furthermore, we update the leakage scheme to include special relativistic 



corrections for the first time to our knowledge (see Appendix B for details). 



Spherical coordinates, (r, 9, cf)) are used in our simulations and the computational domain 
is extended over 50km < r < 30000km and < 6* < 7r/2 and cover ed with 300 (r) x 40 (6^) 
meshe s with the assumption o f axial and equatorial symmetry. As in iMacFadyen fc Woosley 
( I1999I ): iFujimoto et al.l (120061 ) . we adopt an absorptive inner boundary condition at 50 km. 
Together with SR treatment of MHD, this inner boundary, albeit taken to be large, allows us 
to explore the long-term evolution of collapsars. Note that we do not necessarily assume the 
formation of BH inside, but only assume that the central region would not affect the regions 
outside. One interpretation of the position of the inner-boundary could be a surface of the 
standing accretion shock waves produced after bounce. In §5.31 by counting the accreted 
mass in the central objects when the simulations terminate, we will discuss what they could 
be, namely the neutron stars or the BHs. 



The total gravitational potential in Eq. (4) is 



tot 



PW 



POI) 



(7) 



wher e $pw mimics the contribution from strong gravity around the central objects (jPaczynsky fc Wiita 
1980l ). Under the special relativistic modification, this potential has been suggested as use- 
ful to reproduce the dyn a mics outside the last stable orbit in the Schwarzschild metric 
( lAbramowicz et al.l (119961 ): iFukud (120041 )). Thus such treatment, albeit very approximate, 
may not be too bad for our computations (see discussions in §5.31) . To achieve further ac- 
curacy, we need to perform simulations in curved time-space with general relativistic MHD, 
which is beyond the scope of this paper. 



Initial conditions 



As for the initial profiles of the collapsing star, we employ the spherical data set of 
the density, temperature, internal energy, an d electron fraction in m odel 350C with the 
mentioned chemically homogeneous evolution (jWoosley &: Heger 2006). At the the zero age 
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Table 1: Models and Parameters. Model names are labeled by the initial strength of mag- 
netic fields and rotation. Bq is a constant in Eq. fllOp . a is the ratio of the specific angular 
momentum normalized by the one at the last stable orbit in Eq.(IH]). r/|iy| and -Emag/lW^I 
represents the ratio of the rotational energy and the magnetic energy to the absolute value 
of the gravitational energy, respectively. 



main-sequence, the progenitor mass, rotational velocity, and metallicity are 35M0, = 
380km/s, and O.IZ0, respectively. At the presupernova phase, the stellar mass is striped 
to be 28.O7M0 due to the mass-loss and the mass of the central iron core is 2.O2M0. Our 
numerical grid contains inner 8.56M0 of the star. 

Since little is known about the spatial configurations of the rotation and the magnetic 
fields from the stellar evolution calculations assuming spherical symmetry of the stars, we 
add the following rotation and magnetic field profiles in a parametric manner to the core 
mentioned above. 

As for the initial angular momentum of the core, we par ametrize the strength by th e 
angular momentum of the l ast stable orbit (: i isn) following iLee fc Ramirez-RuizI (120061 ): 
Lopez-Camara et all J2009h : IProga et al.l J2003al lbl): IProgal J2OO3I . l2005h as 



J = ajUM{R)), (8) 
where j is the specific angular momentum, M{R) is the spherical mass coordinate, encom- 
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passing the mass inside radius R, and a is a model parameter. In this study, we set 1 < a < 3. 
Note that, ahhough angular momentum is larger than jiso with this range, it does not assure 
the formation of the stable disk because of the existence of the slowly rotating matter in the 
polar regions and relatively large inner boundary of our models. Thus dynamical simulations 
are necessary to specify the criteria of the disk formation. 

As for the initial configuration of the magnetic fields, we assume that the poloidal field 
is nearly uniform and parallel to the rotational axis inside the core and dipolar outside. For 
the purpose, we consider the following effective vector potential, 

Ar = Ae = 0, (9) 

^* = ^^T-^^^i^^' (10) 

where Aj.^g^(p is the vector potential in the r, 6, direction, respectively, r is the radius, tq is 
the radius of the core, and Bq is the model constant. We set tq = 3000 km between the iron 
core and the silicon layers and change parametrically Bq as Bq = 10^, 10^ and 10^*^ G for 
each model. 

We compute 15 models changing the initial angular momentum and the strength of 
magnetic fields by varying the value of a and Bq. Each models are named as BXJY, where 
X indicates the initial poloidal magnetic field (10"^ G), and Y represents the ratio of the 
specific angular momentum to jiso- For example, B9J1.5 represents the model with Bq = 10^ 
and j = l.Sjiso- The model parameters are shown in Table [H It is noted that T/|W^| 
and -Emag/lW^I for the original progenitor of the model 350C is 2 x 10~^ and 1 x 10^^, 
respectively. Thus the model series with Jl.O have almost the same angular momentum with 
the progenitor. Considering the mentioned uncertainties of the progenitor models, we choose 
to explore relatively smaller field strength, which has been less investigated so far. 



4. Results 

Computing 15 models in a longer time stretch than ever among previous coUapsar 
models, we observe a wide variety of the dynamics changing drastically with time. To capture 
the general properties of all the models, we first pay attention to the time evolutions of the 
central mass, the mass of the accretion disk, and the neutrino luminosity in the following. 
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Fig. 1. — Time evolutions of the masses accreting to the central objects (left) and the masses 
of the accretion disks (right). Top 5 lines of red, green, blue, purple, and water correspond to 
the models with the same initial magnetic field of 10^ G but with j = ji^o, l-5jiso, 2jiso, 2.5jiso, 
and 3jiso, respectively. For the case of j = 2jiso (blue), the dashed, solid, and dotted line 
show the variation of the initial magnetic fields {Bq = 10^°, 10^, and 10* G), respectively. 
It can be shown that for rapidly rotating models, the mass accretion rate to the center 
becomes smaller (left) and the accretion disks become heavier (right). To estimate the disk 
mass, we count the mass elements which are nearly in the hydrodynamical equilibrium near 
the equatorial plane. 



4.1. General features 



4.1.1. Central mass and disk mass 



The left panel of Figure [H shows the time evolution of the central mass for some repre- 
sentative models, where the central mass is defined to be the baryonic mass accreted onto 
the central object through the inner boundary. It is shown that the central mass grows larger 
and more quickly for the models with the smaller initial angular momentum. This is simply 
due to the smaller centrifugal forces. In fact, the slowest rotation model of B9J1.0 shows a 
fastest increase, and as a result, the mass of the accretion disk becomes lightest (right panel). 
Such a feature seems close to the so-called dwarf disk, in which matter is swallowed to the 



center with nearly a free-fa ll velocity in the equatorial plane (iLee fc Ramirez-Ruia (120061 ): 
Lopez-Camara et al.l (120091 )). Here it should be noted that very rapid gr owth of the central 



mass has been predic ted to affect the coUapsar ability to produce jets (jJaniuk et al.l 12008 



Janiuk fc Progall2008l ). which we observe actually in our simulations and will be explained in 
the following section. For all the models with Jl.O, on the other hand, we do not observe any 
MHD jets owing to inefficient winding of the magnetic fields. Moreover neutrino luminosi- 
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ties from very thin accretion disks for the models is typically about 2 orders-of-magnitudes 
smaller than for the other models. So, we regard such models as inadequate for the progen- 
itor of collapsars and focus on the models with j > 1.5jiso in the subsequent sections. Here, 
we define the critical angular momentum to be the minimum angular momentum to form 
the stable and thick disk. In our case, the critical value is j = l.Sjiso- 

From the left panel, it can be also seen that the increase of the central mass becomes 
slightly larger for models with larger initial magnetic fields (compare blue lines). This is 
likely due to the angular momentum transport via the magnetic fields, which enhances the 
matter accretion to the center. However the difference of the initial angular momentum is 
more decisive to capture the main dynamical features among the computed models, on which 
we focus in the following. 

4-. 1.2. MHD outflows and neutrino luminosity 

To extract the general features furthermore among the models, we focus on the proper- 
ties of MHD outflows and neutrino luminosities. It is noted that both of them is helpful to 
understand the energy sources for powering the GRBs namely via magnetic and/or neutrino- 
heating mechanisms. More specifically speaking, the information of the neutrino luminosity 
is indispensable to estimate the neutrino energy deposition from the accretion disk. And the 
formations of MHD outflows is also important, because it can evacuate the funnel for the 
secondary jets along the rotational axis, which could be the birthplace of relativistic fireballs. 

Top column of each cell of Table [2] indicates whether the MHD jets are formed (yes or 
no indicated hj Q ot x). TYPE I or II indicates the difference of the formation process of 
the MHD jets, which we will explain in detail from the next section (section l4.2p . 

The quantities of the middle column show the neutrino luminosity (sum of all the 
neutrino species z/g, z/g, and ux) estimated at the epoch when the accretion disks becomes 
stationary(bottom) (e.g., typically ~ 4 sec in Figured]). We find that the neutrino lumi- 
nosities become higher for slower rotation models. This is because the accretion disks can 
attain higher temperatures due to the gravitational compression. It is interesting to note 
that the luminosities tend to become smaller for strongly magnetized models with relatively 
smaller angular momentum (J < J2.5). This is mainly because the gravitational compres- 
sion is hindered by the magnetic forces confined in the disks. More detailed analysis of the 
luminosities is given in section 14.31 
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Model Name 


J1.5 


J2.0 


J2.5 


J3.0 




(TYPE II) 


(TYPE II) 


(TYPE I) 


(TYPE I) 




1.1 X 10^"^ erg/s 


4.5 X 10^^ erg/s 


1.6 X 10^*^ erg/s 


9.0 X 10*^ erg/s 


BIO 


3.4 s 


5.8 s 


7.7 s 


9.3 s 




O (TYPE I) 


O (TYPE I) 


X 


X 




1.6 X 10^^ erg/s 


5.1 X 10^^ erg/s 


1.4 X 10^° erg/s 


2.5 X 10"^^ erg/s 


B9 


9.2 s 


5.3 s 


12 s 


14 s 




X 


X 


X 


X 




1.8 X 10^2 erg/s 


8.5 X 10^^ erg/s 


1.7 X 10^° erg/s 


4.5 X 10"^^ erg/s 


B8 


4.3 s 


6.0 s 


10 s 


12 s 



Table 2: Properties of MHD jets and neutrino luminosities. Contents of each cell is, whether 
the MHD jets are formed (yes or no indicated hj Q ot x) with the different formation 
mechanisms indicated by TYPE I or TYPE II (top), the neutrino luminosity (middle) esti- 
mated at the epoch (bottom) when the accretion disks become almost stationary. While the 
success or failure of the jet-launch depends both on the initial strength of magnetic fields 
and rotation, the neutrino luminosity is shown to be predominantly determined by the initial 
rotation rates. 



4.2. Formation of Magnetically-dominated Jets 

In Table 2, we categorized the launching of the MHD jets in two ways. Here it should be 
noted that we distinguish the collimated outflow as " jets " where their half-opening angle 
is less than 10°. Before discussing the details of each type in section |4.2.2[ we mention the 
amplification of the toroidal magnetic fields and the subsequent formation of the magnetic 
outflows from the accretion disks, which precedes the jet formations. 



4-2.1. Amplification of magnetic fields and the outflows from accretion disks 



Since the initial models investigated here are assumed to have only the poloidal fields 
(section [3]), the key ingredients for amplifying the toroidal fields are the compression of the 
poloidal fields and the efficient wrapping of them via differential rotations. In addition, 
MRI should also play an important role, whose wavelength of the fastest growing mode is 

10"'g/cm^ ' 



1/2 



given by A ~ 5(^^j [t^) (^ "" y"'" j km flBalbus fc Hawlevlll998D . Here, putting the 
typical physical values of the disk, our numerical grids are insufficient to capture MRI at 
earlier phase when the magnetic field is weak, but can handle it in the later phase when the 
magnetic field gets stronger. In this sense, the discussion below, which fails to include the 
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effects of MRI completely, should give a lower bound for the field amplification. Discussions 
about MRI are given with numerical tests in section 15.41 

Figure [2] shows the distributions of density and poloidal fields of model B9J1.5 at 7.99 
s, just before the launch of the MHD outflows from the accretion disk. The density takes 
its maximum value at around 120 km in the equatorial plane (left panel), in which the 
poloidal fields are strong because the higher compression is achieved (right panel). The 
white solid line in the right panel indicates the surface positions of the accretion disk. So 
the amplifications of the poloidal fields occur most efficiently inside the accretion disk. It 
is noted here that the disk is gravitationally stable because the adiabatic index (7) inside 
the disk becomes greater than 4/3 due to the contribution of the non-relativistic nucleon 
(7 = 5/3) photodissociated from the iron nuclei. 

As for the toroidal fields. Figure [3] shows that their amplification rates are highest also 
inside the disk (right), because the degree of the differential rotation is large there (left). 
In previous collaps ar simulations assumin g much larger angular momentum initially (e.g.. 



Proga et al.l (j2003bl ): lFujimoto et all (120061)), it seems to be widely agreed that the differential 



rotation is a primary agent to amplify the toroidal fields. On the other ha nd, ou r resul ts 



show that for long-term evolution of relatively slow rotation models (see also iProgal (120051 )). 
the amplification of the poloidal fields by compression is preconditioned for the amplification 
of the toroidal fields. 

The left panel of Figure H] shows that the toroidal fields are amplified as high as 10^^ G 
inside the disk, where the amplification rate is indeed high (see Figure |3] (right)). The ratio 
of the magnetic to rotational energy at this moment is about ~ 1 0%, which is close to the 



saturation level of the field growth as shown in lShibata et al.l (120061 ). From the right panel, it 



can be seen that the magneto-driven outflows, in which the magnetic pressure dominates over 
the matter pressure, are produced in the vicinity of the accretion disk. Looking carefully, 
the amplification rates are higher near along the equator (Figure 3 (right)) and decrease 
as the distance to the equator gets larger vertically (perpendicular to the equator). We 
find that the vertical gradient of the magnetic pressure near the surface of the disk can 
drive the MHD outflows. It is interesting that the propagation of the outflows is not along 
the rotational axis, but slightly off-axis. This is because the ram pressures just along the 
rotational axis, free from the centrifugal forces, are highest, and thus the magnetic pressures 
cannot overwhelm the ram pressure there. Following these magnetic outflows from the disk, 
the MHD jets are formed in the two ways (namely type I or II), which we will explain from 
the next section. 
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-3e+07 -1.5e+07 1.5e+07 3e+07 

Fig. 2. — Logarithmic contour of density (left) and poloidal magnetic field (right) for B9J1.5 
model, just before launching an outflow at 7.99 s. The white solid line (right) denotes the 
area where the density is equal to 10^^ g/cm^, showing the surface of the accretion disk. 
The strong poloidal fields inside the accretion disk in the vicinity of the equatorial plane 
(right) coincide with the high density regions (left). This means that the poloidal fields 
are amplified mainly by compression. Note that the velocity fields are drawn by the white 
arrows, and the length is normalized by the scale shown in top right edge of the box (here, 
lO^'^cm s~^). The central black circle (50 km in radius) represents the inner boundary of our 
computations. In the following figures, density line, velocity field, and the inner boundary 
are depicted in the same manner. 
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Fig. 3. — Same as Figure 2 but for logarithmic contour of differential rotation {-.XdVt/dX) 
(left) and amplification rate of toroidal magnetic field {:dB^/dt) (right). Comparing with 
the right side of Figure [2], the amplification rates become larger where the poloidal fields and 
the degree of the differential rotation are stronger. 
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Fig. 4. — Same as Figure [3], but for the logarithmic contour of toroidal magnetic fields (left) 
and the inverse of plasma beta (e.g., /3~^g = B ■ B/Svr/p) (right). It is shown that the 
toroidal fields are amplified up to 10^^ G (left), and that the MHD outflows are driven from 
the accretion disk (right). In this figure, we plot the velocity fields only for Vr > lO^cm/s 
for illuminating the outflows, where Vr is the radial velocity. 
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4.2.2. Two types of MHD jet formation 

When the magnetically-dominated outflows mentioned above are so strong enough to 
come out of the central iron cores without shock-stall, we call them as type II jets. Even if 
these prompt outflows stall at first, we find an another way of launching jets (type I), whose 
formation processes are a little bit complicated than for type II. We first explain type II in 
the following. 

Figure [5] shows the evolution of the MHD outflow for model B10J1.5, from its initiation 
near from the inner edge of the accretion disk (top left), propagation along the polar axis (top 
right), till they come out of the iron core (bottom). Note the difference of the length scales in 
each panel. Among the computed models, this model has the strongest initial magnetic fields 
with smallest angular momentum (e.g.. Table 2). The toroidal fields can be much stronger 
than for other models, by the enhanced compression of the poloidal fields inside the disk, as 
mentioned in §4.2.11 As a result, the MHD outflow is so strong that they do not stall, once 
they are launched. In fact, the outflow is shown to be kept magnetically-dominated (inverse 
of the plasma beta greater than 1, right-hand side in Figure [5]) till the shock break-out. As 
the disk becomes more compressed, the magnetic pressure of toroidal fields and its vertical 
gradient inside the accretion disk become larger, which acts to push the outflow vertically 
more strongly. In the early phase of the jets (top left), the sideway expansion of the outflow 
is suppressed by the external pressure Pcxt) which is determined by the ram pressure of a 
freely falling fluid with velocity of t^^cxt Pcxt ~ P^r^ ^ ~ pGM /vcxt, with Text being the width 
of the outflow. Text is determined by the balance equation of P^/Stt = Pcxt, which gives a 
reasonable value. This confinement promotes the outflow to keep progressing vertically. As 
the outflow propagates rather further from the center (top right), the outflow begins to be 
coUimated due to the magnetic hoop stresses, and keep their shape till the shock-breakout 
(Figure [6]). It is interesting to notice that the ram pressure just along the polar axis is so 
large that no outflow is formed there, and that there stays a polar funnel, where the material 
accretes onto the central object. We speculate that the formation of the funnel in such an 
early phase could possibly provide nice environments as a birthplace of fireballs, because the 
neutrino heating from the disk could be sufficiently high at the epoch as will be discussed in 
section 5.2. 

It should be noted here that the above outflow driven by the t oroidal fields is es- 



senti ally same as the 'magnetic tower' which was first introduced by iLynden-Belll ( 11996 



20031) in the context of active ga lactic nucleus, and applied to the c oUapsar environments 



( Uzdensky &: MacFadyenll2007al lbl). However in the analytic models by lUzdensky &: MacFadyen 
( 2007a 3) , the driven mechanism of the tower are assumed to be the winding of the magnetic 



fields threaded in the planar accretion disks with no vertical structures inside. Our simula- 
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tions suggest that the origin comes from the vertical gradient of the twisted toroidal fields 
inside the accretion disk. 
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Fig. 5. — Evolution of MHD jets launched from the accretion disk (Type II jets). In each 
panel, logarithmic contour of entropy (left side) and the inverse of plasma beta (e.g., /3~^g = 
B-B/Stt/p) (right side) are shown at 1.43 s (top left), 1.87 s (top right), and 2.11 s (bottom) 
for model B10J1.5. Note the difference of the length scales among panels. Without shock- 
stall, magnetically-driven outflows (/5~ag ^ 1) come out of the iron core (~ 3000 km in 
radius) . 
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Fig. 6. — Three dimensional plot of density with the magnetic field lines (silver line) for 
B10J1.5 model near at the moment of the shock break-out (2.11 s). Color contour on the 
two-dimensional slice represents the logarithmic density. The outer edge of the sphere colored 
by blue represents the radius of 2 x lO^cm. The outflow is shown to be driven by the so-called 
magnetic tower, i.e., by the toroidal flelds tangled around the rotational axis. Note that the 
field lines outside the bluerish region, seen to be more weakly twisted than inside, come 
mainly from the preshock region. 
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Now we move on to discuss the jets in type I, by taking model B9J1.5 as an example. 
In this case, the magnetic outflow from the accretion disk is not as strong as model B10J1.5, 
and stalls at flrst in the iron core (see butterfly-like regions colored by red in the top left 
panel (right-hand side) in Figure [7]). In the top right panel, very narrow regions near along 
the rotational axis are seen to be produced in which the magnetic pressure dominate over the 
matter pressure (colored by red in the right-hand side). Such regions are formed by turbulent 
inflows of the accreting material from the equator, crossing the butterfly-like regions outside 
the disk, to the polar regions. Such flow-in materials obtain sufficient magnetic amplifications 
when they approach to the rotational axis where the differential rotation is stronger, leading 
to the formations of the MHD outflows along the rotational axis (bottom). 

Figure [S] shows the magnetic pressure of the toroidal flelds and the ram pressure of the 
accreting material near along the rotational axis for the same time epoch as in Figure [71 The 
mentioned turbulent inflows to the polar regions happens at 8.148 s (top right), showing the 
increase of the magnetic pressure in the central region within 100 km. When the toroidal 
flelds are amplifled enough to overwhelm the local ram pressure, pv^ with Vr being the radial 
velocity (bottom), the flow-in materials begin to move outwards (bottom in Figure [T]). In 
fact, at about 70 ms later, the initiated outflows can produce the well-collimated magneto- 
driven explosions, reaching ~ 1000 km along the polar axis (Figure [9]). Although the shape 
of the narrow jets are much like the ones obtained in previous collapsar s imulations with 
much rapid ro tation a nd str o ng magnetic flelds (e.g. , iFuiirnoto et al.l (120061 ) ) or with slower 
rotation (e.g., Proga ( 2005 ): Moscibrodzka fc Proga ( 20091 )). it should be noted that their 



formation processes are different. The obtained outflow here is peculiar as a consequence 
of the long-term evolution of coUapsars, which is produced by the interplay between the 
decreasing ram pressure and the magnetic twisting of the turbulent inflow in the vicinity of 
the polar region. 

Now let us return to Table [2] again. Among the models with outflows, models B10J1.5 
and B10J2.0 make the type II jets, and the other make the type I jets. As mentioned, the 
type II is obtained for models with stronger magnetic flelds with relatively smaller angular 
momentum. If each of these conditions were not satisfled, the MHD outflows would stall at 
one time but with the subsequent revival (type I) or stall forever (x in Table 2). 
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Fig. 7. — Same as Figure 4 but for model B9J1.5 at 8.115 s (top left), 8.148 s (top right), 
and 8.154 s (bottom), showing the moment of the formations of jets in type I. In the top 
right panel, very narrow magneto-driven regions along the rotational axis are produced by 
inflows of the accreting material in the equator to the polar regions. Such flow-in materials 
are shown to start propagating along the rotational axis (bottom). 
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Fig. 8. — The magnetic pressure of toroidal fields vs. the ram pressure for model B9J1.5 
at 8.115 s (top left), 8.148 s (top right), and 8.154 s (bottom) along the rotational axis, 
respectively. Note that the timescales correspond to the ones in Figure [71 Only after the 
magnetic pressure becomes dominant over the ram pressure (bottom), the materials inside 
~ 200 km begin to propagate outwards. 
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Fig. 9. — Logarithmic contour of entropy (left) and inverse of plasma beta (right) for B9J1.5 
model at 8.23 s. The narrow magneto-driven explosions are shown near the rotational axis, 
which is produced by the type I mechanism (see text). High entropy region (~ 20) outside 
the collimated jets (colored by light-blue (left)) is a cocoon, which is produced by a fallback 
of the matter from the shock front. 
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Model 


tjet 


[s] 


Mjet [Mq] 


E^et 


erg] 


Tjet 


r. 

jct,mag 


B10J1.5 


2.66 


1.3 X 10-3 


5.6 X 10^8 


1.0024(0.07c) 


1.043(0.28c) 


B10J2.0 


3.8 





5.1 X 10-^ 


5.5 X 10^^ 


1.00063(0.04c) 


1.0087(0. 13c) 


B10J2.5 


6.33 


9.3 X 10-^ 


6.5 X 10^^ 


1.0021(0.06c) 


1.052(0.31c) 


B10J3.0 


8.8 


9 


1.7 X 10"^ 


2.4 X 10^^ 


1.0077(0.04c) 


1.031(0.24c) 



Table 3: Properties of MHD Jets, tjet is the time when the jets come out of the iron core. 
Mjet and Ej^t is the mass and the explosion energy of jets, respectively. Fjet and Fjet.mag 
represent the Lorentz factors (and the velocity normalized by the speed of light) and their 
maximum values estimated by taking into account the magnetic energies, respectively. Note 
that all of them are estimated at tjet- 



4.2.3. Properties of MHD jets 

Now we proceed to look more in detail to the properties of the MHD jets. Here we focus 
on the models with BIO, because all the model sequence is accompanied with the jets either 
type I or II (Table 2). 

Table [3] shows the mass, explosion energy, and Lorentz factor of jets at the moment 
of shock break-out. Here it is noted that the explosion energy is estimated for the regions 
where the local energy is positive and the radial velocity is positive, indicating that the 



matter is not bound by the gravity (see the definition of ciocai > in Appendix A ). As seen 
in the table, the jet of model B10J1.5 has the largest explosion energy with largest baryon 
loading. This is because the jet is type II as mentioned. Since the jet is launched rather 
earlier (~ 1.9 s) than for type I, there is much material near the rotational axis, which makes 
the baryon-load of jets larger for the model. For type I jets, no systematic dependence of 
the initial angular momentum on the masses and the energies is found. We think that this 
is because the formation of the type I jets occurs by turbulent inflows as already mentioned 
in the previous section. 

The similarity between types I and II is that the jets are at most subrelativistic (0.07c for 
model B10J1.5) with the explosion energy less than lO^^erg. To see the maxim of the Lorentz 
factor in our computations, we boldly assume that the magnetic energy of the fluid is fully 
converted to the kinetic energy, having in mind the dissipative process such as magnetic- 
reconnection (Fjet,mag in Table [3]). Even in this case, the jets become only mildly relativistic. 
While the ordinary GRBs require the highly relativisti c ejecta, we speculate that these mildly 



relativistic ejecta may be favorable for X-ray flashes (jSoderberg et al.ll2006l : iGhisellini et al. 



20071 ). which is a low energy analogue of the GRBs. The propagation of the MHD jets can 
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expel the matter along the polar axis. Such baryon-poor environments may be a favorable 
cite for producing the subsequent jets, which could attain high Lorentz factors pushed by 
the magnetic outflows and/or heated by neutrinos. From the next section, we study the 
properties of neutrino luminosities obtained in our simulations and discuss how the neutrino 
heating, albeit not coupled to the hydrodynamics here for simplicity, could have impacts on 
the jet formations. 
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4.3. Properties of neutrino luminosity 

As seen in Figured] (left), the disk masses for models with the MHD jets (j > l.Sjiso), 
differ several times at most. On the other hand, the neutrino luminosities, which are con- 
tributed mainly from the accretion disk, differ up to 3 orders-of-magnitudes (Table 2). This 
indicates that the local neutrino emissivities also change over several orders-of-magnitudes 
in the vicinity of the disk. In this section, by focusing on the structures of the disk, we 
discuss how the properties of the neutrino luminosities change for the models with different 
initial angular momentum. 

By comparing Figures [TU] with [TTl it can be seen that for smaller initial angular momen- 
tum (Figure [TOl) . the disk becomes more compact (top left) with higher temperatures (top 
right), and that the resulting enhancement of the electron captures lower in the disk (bot- 
tom). The enhanced compression is a primary reason for the higher neutrino luminosities 
for slower rotating models in Table 2. 

The left panel of Figure [12] shows that the disk, whose equatorial size is ~ 200km (top 
left in Figure 10), is very thick to neutrinos. In fact, the opacity inside the disk becomes 
up to 200 for z/g inside the disk, which is one order-of-magnitude larger than that for z/g- 
The right panel shows that the local neutrino emissivities are highly suppressed by exp(— Tj,) 
with being the opacity, which gives the actual neutrino cooling rates in such a thick region 
(dashed lines). The higher suppression for z/g makes its cooling rate almost comparable to 
that of z/g, which is also shown in Figure [13] It is noted that because of this high neutrino 
opacity inside the disk, the characteristic neutrino cooling timescale is typically more than 
four orders-of-magnitudes longer than the advection timescale. Thus the disk of our models 
is advection dominated flow. 

Figure [14] shows the evolution of neutrino luminosities for each neutrino species (z/g (top 
left), z/g (top right), and ux (bottom)). It is shown that for every species, the luminosities 
become larger for models with smaller angular momentum due to the mentioned higher 
compression, and also that the luminosities of z/g and i^g are dominant over the ones of 
vx- For models with J1.5 and J2.0, the luminosities of Ve become almost comparable to 
those of z/g, while the luminosities of z/g are much smaller than those of z/g for more rapidly 
rotating models (compare top left and right panels in Figure 14 noting the difference of 
the vertical scales). On top of the effect of opacity mentioned above, this is because the 
higher compression leads to the production of positrons more abundantly, which promotes 
the production of z/g via the positron captures. 

Figure [TJ] also depicts the effect of the magnetic flelds on the luminosities (compare 
between B8J1.5, B9J1.5, and B10J1.5). For every neutrino species, the luminosities are 
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shown to become smaller for models with larger initial magnetic fields. For stronger mag- 
netized models, the disks expand more due to the magnetic pressure inside, leading to the 
suppression of the contraction. 

Among the computed models, the total neutrino luminosity becomes largest ~ 10^^ erg/ s 
for model with J1.5 series (e.g.. Table 2). It is well known that the energy deposition rate via 
pair neutrino annihilation (i/ + z/ — > e~ + e~^) is proportional to L„ -L,-, wi t h being the neutrino 
and anti-ne utrino luminosities, respectively (jSalmonson fc Wilsorull999l : lAsano fc Fukuyama 



200ll . I2OOOI ). Thus almost the equivalent luminosities of and z/g is advantageous for making 



the deposition rate larger for a given sum of the luminosities. Due to the two factors, the 
models with smallest angular momentum computed here are expected to be most suitable 
for making the fireballs via neutrinos. 
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Fig. 10. — Contour of density [g/cm^](top left), temperature [Me V] (top right), and electron 
fraction (bottom) for model B9J1.5. These panels are for 9.22 s, when the accretion disk 
becomes stationary (see left panel of Figure 1). The white line marks p = 10^^ g/cm^, 
representing the surface of the disk. Comparing with Figure [H] which has larger angular 
momentum initially, it can be seen that the disk here is more compact (top left) with higher 
temperatures (top right), and that the resulting enhancement of the electron captures lowers 
Ye in the disk (bottom). 
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Fig. 11. — Same as Figure fTOl but for model B9J3.0. Note here that the white hue marks 
p = 10^° g/cm^. 
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Equatorial Radius [cm] Equatorial Radius [cm] 



Fig. 12. — Plots of opacity (left), emissivity (solid line, right), and cooling rate (dashed line, 
right) for v^. (red) and i^e(green) along the equatorial plane of model B9J1.5 at 9.22 s. The 
cooling rates are greatly reduced from the local neutrino emissivities for the opaque regions 
(< 200 km). For such regions, the higher suppression for makes its cooling rate almost 
comparable to that of v^- 
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Fig. 13. — Logarithmic contour of the coohng rates for z/g (left), and for z/g (right) for model 
B9J1.5 model at 9.22 s. It is shown that the luminosities of Uf, and are almost comparable. 
Regions colored by black have lower values than the limit of the color legend, which is out of 
our interest here. The white line marks p = 10^^ g/cm^, indicating the surface of the disk. 
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Time [s] 

Fig. 14. — Plots of time evolution of the luminosity for z/g (top left), for (top right), and 
for ux (bottom). For slower models with J1.5 and J2.0, the neutrino luminosities of Ue and 
Ue become almost comparable, which is expected to be advantageous for energetizing the 
jets via neutrino heating. 
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5. Discussion 

Here we shall discuss limitations of our simulations and requirements towards more 
sophisticated numerical modeling of collapsar, such as the effects of inner boundary (section 
15. ip . the effects of neutrino heating (section l5.2p . and the importance of general relativity 
(section [5.3p . The latter two is neglected or treated by a simple approximation in this paper. 
Numerical tests are given in section 15.41 



5.1. Effects of inner boundary 



First of all, we discuss possible drawbacks due to the large inner boundary (50 km) 
taken in this simulation. For the Schwartzschild metric, the marginal stable orbit can exist 
for j > 2rgC where Vg is the gravitational radius. Such orbiting flow can lead to the formation 
of the shock on the equat orial plane, a r id is p r edicted to res u lt in the formation o f stabl e 
and thick disk as shown in iProga et al.l (|2003aJ) ; iProga] (120031 ) ; iLopez-Camara et al.l (120091 ) . 



To capture such a feature needs the inner boundary as small as -Rin > 2rg. However, our 
inner boundary is larger than the value. For example, models B10J1.5 and B9J1.5 have the 
central object with mass of 2.5M0 and 3.5M0 at the end of computation, which corresponds 
to rg = 7.5km, 10.5km, respectively. Our inner boundary corresponds to 6.7rg,4.8rg for 
each model. Therefore the critical angular momentum to form the stable accretion disk 
discussed in this paper could be affected by the position of the inner boundary. Moreover 
the large inner boundary, which excises the inner edge of the accretion disk, should lead to 
the underestimation of the neutrino luminosity and the resulting neutrino heating, which we 
will discuss in the next section. To clarify those points, we are now planning to implement 
more compact inner boundary in the long-term evolution, which is computationally very 
demanding, thus sequel of this paper. 



5.2. Importance of Neutrino Heating 

We try to estimate the effects of the neutrino energy deposition via neutrino pair anni- 
hilation (i/ + z/ — > e~ + e"*") from the accretion disk to the polar funnel. The neutrino heating 
is important, however, not included in the simulations here for simplicity. 

In the following, we present an order-of-magnitude estimation of the heating rate. To do 



so, we derive the heating rate with the special relativistic corrections (see Appendix C). For 



simpli city, we take the so-called optically thin limit in the accretion disk as in lAsano Sz Fukuyama 



( 12001! ) and consider the neutrino heating only along the rotatinal axis. It is noted that ex- 
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pect for the polar funnel, to make neutrino-heated outflows is hopeless due to the baryon 
contamination. By comparing timescales such as heating and advection, we discuss how 
important the heating could be. 

Figure [15] shows various timescales. Trei = p(? jq^ characterises the timescale for matter 
to become relativistic by the neutrino heating, in which p and represent the local density 
and the heating rate, Trd = p(? jq^ is the timescale when the motion of the matter becomes 
relativistic due to the energy deposition. Tint = ^mtlq^ is the timescale when the neutrino 
heating is comparable to the internal energy (cint) and thus affects the dynamics, Thyd = 
X/I^rl indicates the hydrodynamical timescale with X and Vr being the length and radial 
velocity along the polar axis, rught = Xjc is the light crossing timescale. As shown. Tint 
gets shorter to be several milliseconds near around 100 km, where the timescales become 
most close to Thyd- This means that rather in the vicinity of the center (< 200 km along the 
rotational axis), the neutrino heating have potential importance to affect the hydrodynamics. 
However, it should be noted that it does not directly mean that the heated matter could 
become relativistic. In fact, t^c\ is at least two orders-of-magnitudes larger than rught- 

For making the outflows relativistic, one possible way is to decrease Trei by lowering 
the density of the funnel regions. Such low density regions would be formed if we continue 
to follow the dynamics of coUapsar in more longer term. Unfortunately however, the nu- 
merical difficulty of treating such force-free fields prevents us from doin g so. The ii u meric al 



code specially developed to solve the force-free fields is required (e.g., iMcKinneyl (l2006l )). 
which is major undertaking. Another way is to increase the heating rate q^ . If the neu- 
trino heating not just from the equatorial plane here but also from the who le accretion disk 
could be included such as by the ray-tracing calculation (IBirkl et al.l 120071 ) . the deposition 
rates would become larger due to the ge ometrical effects. D issipative processes such as 
magnetic-reconnection/ Joule heating (e.g., iProga et al.l (j2003al )) inside the disk would rise 
the temperature, which could be also good for achieving the higher luminosity. Further- 
more, general relativis t ic effe c ts would increase the depos i tion r ate in the vicinity of the BH 
dSalmonson fc Wilson! JlQQoh : lAsano fc Fukuvamal J200ll . boooh ). which is also remained to 
be studied. 



5.3. General relativistic effects 



As mentioned earlier, we have employed the Paczynsky-Wiita potential to mimic the 
Schwarzschild metric. With the special relativistic modification, this artificial potential is 
known to be able to approximate the general rel ativistic (GR) motion well for the r egion s r > 
"ivg (with Tg being the Schwarzshild radius, e.g.. lAbramowicz et al.l (ll996l ). lFukud (120041 )). to 
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which we focused on in this study. Howeve r needless to say, wh at is needed for coUaps ars is 



GR simulations a round the Kerr BH (e.g. , IShibata et all (120071 ): iMontero et al.l (120081 ) and 



see discussions in IWoosley fc Hegerl (120061 ) ). 



The maximum mass of neutron sta rs, depending on nuclear equations of state, i s esti- 
mated to be less than typically 3.OM0 ( iLattimer fc Prakaslj|200ll : IZhang et al.ll2008l ). Fol- 
lowing this criteria, the central objects of models with smaller initial angular momenta of 
Jl.O and J1.5 may collapse into the BHs near within 3 s (Figured] (left)). In such cases, GR 
effects very close to the inner edge of the accretion disk should be important and their im- 
pacts on the MHD outflows and the neutrino heating should be investigated. It is naturally 
expected that strong gravity due to GR effects will lead to not only the efficient gravitational- 
wave emission, but also the enhanced neutrino emission due to the cor apression. As r ecentl y 



studied extensively in the context of core-collapse supernovae (e.g., iKotake et al.l (120091 ): 
Kawagoe et al.l (120091 ): lOttI (120091 ) and references therein), gravitational-wave and neutrino 
signatures also from coUapsars should give us a new observational window to probe the cen- 
tral engine. This paper is a prelude before our forthcoming work to clarify those aspects, 
which will be presented elsewhere soon. 



5.4. Numerical Tests 

Figures [16] and [17] show a convergence of our numerical results for different grid res- 
olutions. Figure [16] shows an agreement of the neutrino luminosity. This means that the 
evolution of the accretion disk, whose temperature and density profiles determine the neu- 
trino luminosity, is numerically converged for the different resolution we tested. On the other 
hand. Figure [T7| shows that the magnetic energies before ~ 3 sec are rather sensitive to the 
numerical resolution, although the overall trends are similar. We suspect that the discrep- 
ancy in the earlier phase should come mainly from the effects of MRI, which are difficult to 
be captured by our numerical resolution as already mentioned in 14.2. 1[ 

Then we discuss the validity of the equatorial symmetry assumed in this study. Figure 
[T8] shows that the obtained luminosities agree well with each other, however the magnetic 
energy differs up to two times (Figure [T9]) . This could be also due to MRI. As mentioned, our 
numerical resolutions can treat MRI marginally in the sense that it can follow MRI in the late 
phase when the field strength becomes stronger, while it cannot resolve MRI sufficiently in 
the early phase when the field strength is weaker. In this sense, the obtained results should 
give a lower bound for the criterion of the field amplification and the jet formation. To 
capture MRI fully is unquestionably important for collapsar simulatioi is, but possibly needs 



the prescription such as by the adaptive mesh refinement scheme (e.g.. lZhang fc MacFadyen 



( 120091 )) , which we pose as a future task. 



6. Summary 

In hght of the coUapsar model of gamma-ray bursts (GRBs), we presented our numeri- 
cal results of two-dimensional magnetohydrodynamic (MHD) simulations of the collapse of 
rotating massive stars. Pushed by recent evolution calculations of GRB progenitors, we fo- 
cused on lower angular momentum of the central core than previously assumed. As for the 
initial magnetic field strength, we chose to explore relatively smaller field strength (< lO^'^ 
G), which has been less investigated so far. By performing special relativistic simulations 
including both realistic equation of state and neutrino cooling, we followed a long-term evo- 
lution of the slowly rotating collapsars up to ~ 10 s. Then we studied how the formation 
of MHD jets and the properties of accretion disks could change with the initial angular mo- 
mentum and the initial magnetic field strength. The obtained results can be summarized as 
follows. 

1. Our numerical results show that for the GRB progenitors to function as collapsars, 
there is a critical initial angular momentum: jcrit — l.Sjiso with jiso being the angular 
momentum of the last stable orbit, below which matter is quickly swallowed to the 
central objects, no accretion disks and no MHD outflows are formed. 

2. When larger than the criteria, we find that smaller initial angular momentum leads 
to more compact accretion disk due to compression. It seemed widely to be agreed 
in previous collapsar simulations that the differential rotation is a primary agent to 
amplify the toroidal fields and the resulting MHD outflows. On the other hand, our 
results show that for relatively slow rotation models, the amplification of the poloidal 
fields by compression is preconditioned for the amplification of the toroidal fields and 
the MHD outflows. 

3. Among the computed models, we find the launch of the MHD jets in the following two 
ways. For models with stronger initial magnetic fields {Bq > 10^° G), the gradient of 
the magnetic pressure perpendicular to the equatorial plane inside the accretion disk 
can drive the MHD outflow. This outflow makes the strong magnetic explosions like a 
'magnetic tower', which we called as type II. For models with weaker initial magnetic 
fields, the magnetic tower stalls first and the subsequent MHD outflow is produced by 
the accreting material from the equator to the polar region. Type I jet is found to be 
produced when such flow-in materials obtain sufficient magnetic amplifications, due to 
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the strong differential rotation near the rotational axis, to overwhelm the ram pressure 
of the accreting material. 

4. Regardless of type I or II, the jets can attain only mildly relativistic speeds with the 
explosion energy less than lO^^erg. Such events could possibly be related to the X-ray 
flashes. After the passage of the MHD jets, the baryon-poor environments will be left 
behind. Such polar funnels could provide a favorable cite for the subsequent jets, which 
could attain high Lorentz factors pushed by the magnetic outflows and/or heated by 
neutrinos. 

5. To obtain stronger neutrino energy deposition in the polar funnel heated from the 
accretion disk, we flnd that the smaller initial angular momentum is more appropri- 
ate. This is because the gravitational compression makes the temperature of the disk 
higher. When the accretion disk settles to the quasi-stationary state in their late time 
evolution, the maximum neutrino luminosity is found to reach ~ 10^^ erg/s. Due to 
the high neutrino opacity inside the disk, the luminosities of z/g and z/g is found to be- 
come almost comparable, which is advantageous for making the deposition rate larger. 
Based on an order-of-magnitudes estimation of the energy deposition, it is suggested 
that the neutrino heating could be as efficient as the magnetic mechanism to energe- 
tize the outflow. Among the computed models here, the model with the initial angular 
momentum of jcrit ~ l-5jiso and with initial magnetic fleld strength of Bq > 10^'^ G, 
provides the most plausible condition for collapsar, because such models are appropri- 
ate not only for producing the MHD outflows quickly by the magnetic towers, but also 
for obtaining the stronger neutrino heating in the evacuated polar funnel. 
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Appendix A. Definition of energy 



According to iTakiwaki et al.l (120091 ) . we define special relativistic description of local 
energies as follows, 



Ckin 


= pW{W-l), 


(Al) 


Cint 


= eW^+p{W^-l), 


(A2) 




B^W 1 \ 
471 V 21^2 J ' 


(A3) 


^mag 


i 


(A4) 


Clocal 


= Ckin + Cint + Cmag + $tot 


(A5) 



where Ckin is the total kinetic energy, Crot is the rotational energy, Cint is the internal energy, 
cb; is the magnetic energy of i th component, Cmag is the total magnetic energy, and eiocai is 
the total local energy. We use these values to compute the explosion energy of jets in §11 



Appendix B. Special relativistic treatments for neutrino cooling 



Following the procedure in iMihalas fc MihalasI (119841 ). we derive the neutrino cooling 
rates including special relativistic corrections. Neutrino emissivity is often calculated in the 
rest frame of relativistic fluids because radiation is usually isotropic there. When converting 
quantities in the rest frame to the laboratory (lab) frame in which the fluid variables are 
defined and evolved by the hydrodynamic equations, we need to take into account corrections 
from the Lorentz transformations between the two frames. 



The cooling rate in the lab frame is calculated as follows. 



ri{Q, e)dedQ, 



(Bl) 



by summing up the neutrino emissivity of ?7(fi, e)(erg/cm^/s/str) of a neutrino energy of e 
both over a solid angle dfl and over de in the lab frame. To satisfy the number conservation 
under the Lorentz transformation, [j] / e)dtdV dedVL is a Lorentz invaria nt. It is noted that 
the Lorentz transformation can be written in the following simple form (IMihalas fc Mihalas 
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1984h . 



dtdV 
de 

dn 

e 



dtodVo 
—deo 



-0 



dQn 



7(1+ /5/io) Co 
6o/[7(l-/5/i)], 



(B2) 
(B3) 

(B4) 
(B5) 



where /3 = \v\/c,'y = (1 — , and jj = v ■ n/\v\ where n is the unit vector of the 

direction of the emitted neutrino. Here the subscript denotes the variables which are 
measured in the rest frame. Then, we obtain the transformation of the emissivity, 



'n{n,e) = ^?7o(/Uo,eo). 
^0 

When the emission is isotropic in the rest frame, namely as 

?7o(Aio,eo) = Co(eo), 



(B6) 



(B7) 



where Co(eo)[erg/cni'^/s] is the direction-averaged emissivity, the cooling rate in the lab frame 
can be given as. 



ri{Q, e)dedQ 



—VoifJ'O, eo)deodQo 

^0 



7(1 +/3/io) dno 



Co(eo)c?eo 



4vr7 J Co(eo)c?eo- 



7% 



:b81 



This equation means that the cooling rate becomes 7 times larger than the one in the rest 
frame with the special relativistic effect. It is noted that this relation holds when 77o(/io;fo) 
is isotropic in the rest frame, such as inside the accretion disk investigated here, which is 
very opaque to neutrinos and thus the radiation there is well approximated to be isotropic 
(see discussions in section 4.3). 

Now we move on to consider the corrections to the opacity. Again from the number con- 
servation with the Lorentz transformation of the absorption coefficient ^) (erg/cm^/s/str). 
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{x/^)dtdVdedQ is a Lorentz invariant. Therefore, with the same approach in the case of the 
coohng, the absorption coefficient in the lab frame becomes, 



Then the opacity can be calculated as. 



Xo(Ato, eo) 



7(1 - Pl^) Xo(/Uo, eo). 



(B9) 



= 7(1 -/5/i)xo(/io,eo)c?s 
= 7(1- Pfi) dro, 



(BIO) 



where ds is the spatial distance in the lab frame. Here (7(1 — /5/i)) reflects the special 
relativistic effect. 



With the local cooling rate (Eq. fIBSP ) and the opacity (Eq. flBlOp ). we evaluated the 
effective cooling rate appeared in Eq. as 



q exp 



7% exp 



dr 



7(1- /9/i) dro 



(Bll) 



where the integration is performed along the radial direction for simplicity. 



Finally, the reaction rate is also modified in the lab frame by the time contraction. The 
time evolution of Ye is then given by 



dYe 



Todto 
To'jdt, 



(B12) 



where Fq represents the local neutrino reaction rates. Apparently, this modification could 
be important in the regions where the motion becomes (special)relativistic such as the inner 
most region of the disk or inside of the relativistic outfiow. 



Appendix C. Special relativistic modification for neutrino pair annihilation 

Employing the same procedure as in the previous section, we here derive the heating 
rate via neutrino pair annihilation [u + 1? ^ e~ + e~^) with special relativistic corrections. 
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By putting the Lore ntz transformations described in Eq.flBSD into the hea ting rate de- 
rived by previous studies fISalmonson fc Wilsonlll999l : lAsano fc Fukuyamall200ll ). the heating 
rate in the lab frame is given by 



dtdV 



2cKGp J d6„d(f)„d6pd(j)pde^depelel{eu + ep)fu{r,p^,)fp{r,pp) 
X [1 — sin 6,y sin 6p cos {ip,y — ifp) — cos 61, cos 6p]^ sin 61, sin 6p 
2cKGl J de^d^^ddpd(t)u 

+SRt{n,)SRUnp)N,,o{r, fi.)5p,o(r, Qp)] 
X [1 — sin 9i, sin 9p cos {(pi, — ipp) — cos cos 9p]^ sin sin 9p, 



where 



SR,{r,Q,) 



l/[7,(l-/i,/3,)], 
\v,,\ 



1 



Pu 

\Pv\ 



(sin COS 0^, sin sin 01,, cos6'iy) 



(CI) 



(C2) 
(C3) 

(C4) 

(C5) 

(C6) 

(C7) 
(C8) 



and r^^v,^, p^, fy denotes the position of the neutrino source, the velocity of fluid at the 
neutrino source, the momentum vector of neutrino, and the distribution function of neutrino. 
Definitions are same for the anti-electron type neutrino by changing the notation v io v 
from Eq. (lC2p to Eq. llCSp . Subscript again denotes variables which are measured in the 
rest frame of the neutrino source. Here the neutrino source indicates the accretion disk. The 
neutrino number flux along its path to the target is assumed to be conserved for simplicity 
as fy{r,p^) = fviryfl.p^o). In Eq. llCll) . the factor SR^ reflects the special relativistic 
modification to the heating rate. 



The neutrino distribution function inside the accretion is well approximated by the one 
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in the (3 equilibrium as, 

1 driQ 



(hc)'^ deodilodtodVo 
1 1 



where Ty^, rj^^ is the temperature and degeneracy parameter of neutrino in the rest frame, 
respectively. We approximate T^^ to be equal to the temperature of fluid T{r^). Then, S 
and in Eq. (CI) can be expressed by the Fermi integrals JF^ as 



k 

X 



^k{y) = / — —rdx, (CIO) 

J exp(x-i/) + l 



5.,o(r,fi) = ffl^^4(ry,,o), (Cll) 
iV.,o(r,fi) = ^^^.'f;}^' T^M- (C12) 
With these modifications, we calculated the heating rate in section I5.2[ 
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Fig. 15. — Plots of various timescale for matter ; the timescale to get relativistic by the 
neutrino heating Trei (sohd), the timescale to affect dynamics with neutrino heating Tint 
(dashed), the timescale of dynamical motion Thyd (dotted), and the timescale of light crossing 
Tiight (small dotted), at 9.22 s when the accretion disk reaches to the stationary state with 
the funnel regions along the pole. 
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Fig. 16. — Time evolution of neutrino luminosities of (solid line), (dashed line), and ux 
(dotted line) for model B9J1.5. Models labeled by "a" and "b" have doubled mesh points in 
the radial and azimuthal directions than the canonical one of 300 (r) x 40 (^) mesh points. 
Each luminosity coincides well with each other, showing the numerical convergence of the 
obtained results. 
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Fig. 17. — Same as Figure fT6] but for the time evolution of magnetic energy. 
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Fig. 18. — Same as Figure [16] but with and without equatorial symmetry (indicated by 
"sym" and "no sym", respectively). 
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Fig. 19. — Same as Figure fTSl but for the time evolution of magnetic energy. 



